function [betas_slope,betas_yld,R2s_slope,R2s_yld] = yield_regs_nom(gesol, X_Nxz, Y_grid, stationary_prob_NxYz, FLAG_RUN_PARALLEL)
% regresses yld(t) on X(t) for maturities=1,2,...,5 years
% regresses yld(m,t)-yld(1) on X(t) for m=2,...,5 years

    if nargin<5; FLAG_RUN_PARALLEL = false; end

    betas_slope = zeros(1,4);
    R2s_slope = zeros(1,4);
    betas_yld = zeros(1,5);
    R2s_yld = zeros(1,5);

    [NN, Nx, NY, Nz] = size(stationary_prob_NxYz);
    Y_NxYz = permute(repmat(Y_grid(:),[1 NN Nx Nz]),[2 3 1 4]);

    X_NxYz = permute(repmat(X_Nxz,[1 1 1 NY]),[1 2 4 3]);

    % yield regressions
    if FLAG_RUN_PARALLEL
        
        P_A_list = zeros(NN,Nx,Nz,5);
        P_C_a_list = zeros(5,1);
        
        for i = 1:5
            idx = find(gesol.nom_bonds.maturities==12*i);
            P_A_list(:,:,:,i) = gesol.nom_bonds.P_A(:,:,:,idx);
            P_C_a_list(i) = gesol.nom_bonds.P_C_a(idx);
        end
        
        parfor i = 1:5
            yld_NxYz = permute(repmat(-log(P_A_list(:,:,:,i))/i,[1 1 1 NY]),[1 2 4 3]) ...
                + (P_C_a_list(i)/i)*Y_NxYz;
            [betas_yld(i),R2s_yld(i)] = reg_Y_Nz_on_X_Nz(yld_NxYz,X_NxYz,stationary_prob_NxYz);
        end
        
    else
        
        for i = 1:5
            idx = find(gesol.nom_bonds.maturities==12*i);
            yld_NxYz = permute(repmat(-log(gesol.nom_bonds.P_A(:,:,:,idx))/i,[1 1 1 NY]),[1 2 4 3]) ...
                + (gesol.nom_bonds.P_C_a(idx)/i)*Y_NxYz;
            [betas_yld(i),R2s_yld(i)] = reg_Y_Nz_on_X_Nz(yld_NxYz,X_NxYz,stationary_prob_NxYz);
        end
        
    end

    % slope regressions
    idx = find(gesol.nom_bonds.maturities==12);
    yld1_NxYz = permute(repmat(-log(gesol.nom_bonds.P_A(:,:,:,idx)),[1 1 1 NY]),[1 2 4 3]) ...
        + gesol.nom_bonds.P_C_a(idx)*Y_NxYz;
    
    if FLAG_RUN_PARALLEL
        
        P_A_list = zeros(NN,Nx,Nz,4);
        P_C_a_list = zeros(4,1);
        
        for i = 2:5
            idx = find(gesol.nom_bonds.maturities==12*i);
            P_A_list(:,:,:,i-1) = gesol.nom_bonds.P_A(:,:,:,idx);
            P_C_a_list(i-1) = gesol.nom_bonds.P_C_a(idx);
        end
        
        parfor i = 2:5
            yld_NxYz = permute(repmat(-log(P_A_list(:,:,:,i-1))/i,[1 1 1 NY]),[1 2 4 3]) ...
                + (P_C_a_list(i-1)/i)*Y_NxYz;
            [betas_slope(i-1),R2s_slope(i-1)] = reg_Y_Nz_on_X_Nz(yld_NxYz-yld1_NxYz,X_NxYz,stationary_prob_NxYz);
        end
        
    else
    
        for i = 2:5
            idx = find(gesol.nom_bonds.maturities==12*i);
            yld_NxYz = permute(repmat(-log(gesol.nom_bonds.P_A(:,:,:,idx))/i,[1 1 1 NY]),[1 2 4 3]) ...
                + (gesol.nom_bonds.P_C_a(idx)/i)*Y_NxYz;
            [betas_slope(i-1),R2s_slope(i-1)] = reg_Y_Nz_on_X_Nz(yld_NxYz-yld1_NxYz,X_NxYz,stationary_prob_NxYz);
        end
    
    end

end